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ABSTRACT 

We study the properties of massive galaxies at an average redshift of z ~ 0.34 through stacking more 
than 42000 images of Luminous Red Galaxies from the Sloan Digital Sky Survey. This is the largest 
dataset ever used for such an analysis and it allows us to explore the outskirts of massive red galaxies 
at unprecedented physical scales. Our image stacks extend farther than 400 kpc, where the r-band 
profile surface brightness reaches 30 mag arcsec"^. This analysis confirms that the stellar bodies of 
luminous red galaxies follow a simple Sersic profile out to 100 kpc. At larger radii the profiles deviate 
from the best-fit Sersic models and exhibit extra light in the g, r, i and z-band stacks. This excess 
light can probably be attributed to unresolved intragroup or intracluster light or a change in the light 
profile itself. We further show that standard analyses of SDSS-depth images typically miss 20% of 
the total stellar light and underestimate the size of LRGs by 10% compared to our best fit r-band 
Sersic model of n=5.5 and re=13.1 kpc. If the excess light at r>100 kpc is considered to be part of 
the galaxy, the best fit r-band Sersic parameters are n=5.8 and re=13.6 kpc. In addition we study the 
radially dependent stack ellipticity and find an increase with radius from e=0.25 at r=10 kpc to e=0.3 
at r=100 kpc. This provides support that the stellar light that we trace out to at least 100 kpc is 
physically associated with the galaxies themselves and may confirm that the halos of individual LRGs 
have higher ellipticities than their central parts. Lastly we show that the broadband color gradients 
of the stacked images are flat beyond roughly 40 kpc, suggesting that the stellar populations do not 
vary significantly with radius in the outer parts of massive ellipticals. 

Subject headings: galaxies: interactions — galaxies: evolution — galaxies: elliptical — galaxies: 
structure 
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1. INTRODUCTION 

Elliptical galaxies dominate the galaxy mass func- 
tion for M^, > IO-'^^Mq and are the brightest ex- 
tended objects in the nearby Universe. Extensive ef 
fort has been made to study nearby ellipticals, re- 
lying on observations of their morphology, kinemat- 
ics and stel lar populations fe.g.. Faber fc Jackso n! 1 19761 
IllingworthI [T977.: Boroson et al. .1983.: .Dressier et al 



19871: IKormendv fc Diorgovskil Il989l: IPeletieri 119891 
Worthev et al.lll992D . Most studies have focused on the 



bright centers of these galaxies, out to 1-2 effective radii, 
as detection of their faint outskirts has posed a chal- 
lenge both observationally and analytically. Neverthe- 
less, there exists great interest in correctly analyzing 
the full physical extent of massive red galaxies, partic- 
ularly as they may bu i ld up inside-out through mergers 
(iLoeb fc Peeblesll2003l : iNaab et al.ll2007l: iBezanson et al.l 
l2009f ). An accurate measurement of the light profile 
shape of low redshift galaxies is crucial for the interpre- 
tation of high redshift galaxy observations. 

Direct observations of the outskirts of individual mas- 
sive galaxies are hard to perform given the extremely 
faint surface brightness level that is reached at large 
radii. The observed light in such data is highly in- 
fluenced by nearby obj ects, flat flelding errors and the 
wings of the PSF re.g..lMihos et al.ll2005l: Ide Jongll2008l: 
iTa et al.ll20"09tlvan Dokkum et al.ll2009D . A recent study 



Kormendv et all (|2009f ) utilized a compilation of data 
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from several sources including observations of 28 Virgo 
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ellipticals. The authors detected light outside of five ef- 
fective radii in only two galaxies, M87 and M49, which 
are the most luminous Virgo galaxies. There have also 
been many studies of extended light distributions in mas- 
sive clusters, and in those environments the light at 
large distances from the central galaxy is usually at- 
tributed to an intracluster stellar population distinct 
from the central galaxy (Intra-Cluster Light, or ICL). 
The ICL has been successfully t raced to distances of 
>500 kpc in several studies (e.g., [Gonzalez et al.l [20001 : 
IKrick fc Bermteh]l2007l : [Mihos et al.ll2005l) ! 

As an alternative to deep observations, stacking a large 
number of images of similar galaxies could greatly im- 
prove the reached overall depth at the cost of losing sys- 
tem specific information. This technique has been used 
by several author s in studies of both spiral and ellipti- 
cal g al axies (e.g., iZibetti et al.ll200^ Ivan Dokkum et all 
[20101 ). IZibetti et al.l (|2005l hereafter Z05) stacked 683 
images of clusters in SDSS, and reported that the ICL 
accounts for 10% of the total light in galaxy clusters. 

In this study we stack an unprecedentedly large sample 
of LRG images, in a similar way as Z05, to explore the 
faint outskirts of massive red galaxies. LRGs are thought 
to be mostly group centrals that live in lower mass halos 
than the objects studied in Z05. By stacking the LRG 
data we hope to shed new light on the properties of these 
objects at very large radii and to better constrain the size 
and total luminosity of elliptical galaxies at z ~ 0.34. 

2. LRG IMAGE STACKING 
2.1. Sample selection 
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Figure 1. The redshift distribution of tiie entire LRG candidate 
sample (black contour) overlaid with the cuts applied in this study 
(red area). The mean value of the selected sample is z ~0.34. The 
resultant r-magnitude distribution is also shown, along with the 
overall LRG r-magnitudc distribution, in the top-right corner. 



We selected galaxy images for this study from the 
Sloan Digital Sky Survey (SDSS, Abazajian ct al. 2009) 
including all objects classified as Luminous Red Galax- 
ies (LRG) that have a spectroscopic redshift measure- 
ment. LRGs are intrinsically red and luminous objects 
that were identified as such from their central surface 
brightness and locat ion on a rotated color- color diagram 
(for full details see lEisenstein et al.l 120011) . This selec- 
tion is aimed at finding the most luminous red galaxies 
in the nearby Universe [L > 3L*) out to a redshift of 
z — 0.5. Being some of the most massive galaxies in 
SDSS, LRGs occupy the high end of the stellar mass 
spectrum between lO^^Mg and a few times lO^^Afg. 
Roughly 90% of all LRGs are central halo galaxies and 
they mainly reside i n groups with a typical halo mass of a 



few times 10"M„ (jWake et al.ll2008l : iZheng et al.|[200l 



IReid fc SDereel[200a . 

The seventh data release of SDSS (DR7) includes 
188366 spectroscopic LRGs over a wide range of redshifts 
and apparent magnitudes. The redshift distribution has 
two main populations peaking at z ^-^0.05 and z ~0.34. 
The lower redshift LRG candidates are predominantly 
a contamination sub-sample of fainter, lower mass red 
galaxies. In order to assure that our sample is indeed 
composed of high mass, luminous red galaxies we selected 
objects with a narrow distribution of redshifts around the 
high-z peak. This selection also ensures that galaxies in 
this distance and brightness ranges do not suffer from 
significant size and mass evolution within the sample. 

The main sample used in this study comprises 55650 
galaxies in a redshift range 0.28<z<0.40 with mean red- 
shift <2:>=0.34, resulting in an apparent r-magnitude 
range of 18.5±0.4 (figure [J). In fact, 99% of the 
samp le falls within t he flu x limit of cut I, as defined 
by Eisens tein et al.l ()2001l) . making it approximately 
volume-limited. From this master list we excluded 12750 




Figure 2. Image preparation for stacking: thumbnails of size 
200"x200" pixels were cut around each LRG while all other ob- 
jects in the field were masked out. The stack in the top-left corner 
was made using more than 42550 LRG images and can be traced 
to a radius greater than 100 kpc. 

galaxies (23%) where more than 75% of the central 
5" X 5" had to be masked out due to close proximity to 
another object (masking details in subsection I2.2|) . In 
addition, we excluded 293 (<1%) of the LRGs because 
of varying sky levels in the frame caused by close proxim- 
ity to a bright star in or just outside of the field. Finally, 
we excluded 28 (<0.1%) images of galaxies with apparent 
r-magnitude outside of the selected range (details in sub- 
section [2]3]) . The final sample consists of 42579 galaxies. 

2.2. Preparing the images for stacking 

We acquired imaging data for the fields containing the 
selected galaxies from the SDSS archive in all five bands: 
u, g, r, i and z, corresponding to central wavelengths of 
355. Inm, 468. 6nm, 616. 5nm, 748. Inm and 893. Inm, re- 
spectively. For each selected object we cut out a square 
region of roughly 200" x200" (950 kpc at z = 0.34), cen- 
tered on the galaxy, from the five ugriz field images. We 
then shifted the resulting thumbnails using cubic con- 
volution interpolation to center the main galaxy on the 
central pixel. Parts of the resulting thumbnails which 
extended beyond the SDSS stripe edge where given zero 
weight in the stacked images. 

In order to detect and mask out any foreground and 
background objects we created a masking template by 
combining the thumbnails of three of the optical bands 
(r, i and z). This increased the signal-to- noise ratio of the 
template by a factor of roughly compared to the in- 
dividual frames, thus enabling us to unv eil more sources 
in th e field. We then ran SExtractor (iBertin fc Arnout^ 
I1996D on the combined image and detected all the ob- 
jects in the frame. We set the detection threshold to 
a value of 1.4 times the standard deviation above the 
background RMS level and used AUTO photometry to 
extract Kron radii for the detected objects. Finally we 
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Figure 3. Radial ligiit profile of tlie i-band stack (black solid 
line) overlaid with light profiles of the synthetic PSF stack (green 
dotted line), the bright star PSF stack (red dashed line) and the 
combined PSF (solid blue line). Light scatter due to the PSF is 
most significant in i-band images. 



constructed elliptical object masks by growing the semi- 
major and semi-minor axes as determined by Sextractor 
by a factor of 2.5. We found the optimal values for both 
the detection threshold and the mask growth factor from 
trial and error and verified that no additional unmasked 
objects remained in a visual inspection of the images. 
The resulting masks cover on average 18% of the total 
area of each thumbnail. We note that galaxies below 
the detection threshold of SDSS are not excluded by our 
procedure. We return to this in subsection 12.41 

2.3. Magnitude-bin stacks 

After creating object masks using the combined g-(-r-|-i 
images we continued on to preparing the individual band 
images for stacking. We did so by applying the mas- 
ter mask to the images in each of the bands and care- 
fully subtracting both the soft bias, as introduced by the 
SDSS pipeline, and the remaining sky from the images. 
We determined the background sky level by temporarily 
masking out the central galaxy and fitting a Gauss curve 
to the pixel distribution of each masked thumbnail, as- 
suming the residual sky pixel values are distributed like 
Poisson noise. This resulted in a set of five thumbnails 
per field with the central LRG being the only light source 
in each sky-subtracted frame. We compared the back- 
ground level values we found to those obtained from the 
image headers, where those were available, and found 
excellent agreement. 

One of the challenging aspects of stacking images with 
any range of object brightnesses lies in the normalization 
of individual frames in respect to one another. To min- 
imize normalization errors and optimize the S/N ratio 
we did not normalize each individual image but rather 
divided the sample into thirteen bins based on the ap- 
parent magnitude of the galaxies as determined by the 
SDSS pipeline. By doing so we ensured that the galax- 



ies in each bin are "pre-normalized" . The bins span a 
brightness range of 17 to 19.6 r-magnitudes at intervals 
of 0.2 magnitudes, corresponding to a flux difference of 
up to only 20% within any given bin. These magnitude 
cuts were chosen such that only bins that had at least 
100 galaxies in them will be processed. We note that we 
did not rescale the images to a common physical scale 
prior to stacking them. The ±10% variation in physical 
scale over the redshift range of the sample is likely small 
compared to the variation in effective radii between the 
individual galaxies. 

In addition to stacking masked, LRG centered images 
within each magnitude bin we summed up their corre- 
sponding object masks to serve as weight maps for the 
stacks. We then divided the summed images in each 
of the five bands by their respective weight mask, thus 
creating an averaged exposure-corrected stack. Since no 
normalization was required for images within the same 
magnitude bin, the noise characteristics were essentially 
improved by the square-root of the number of images in 
that bin. 

2.4. Random Stacks 

Correct sky subtraction is crucial for properly analyz- 
ing the faint outskirts of the galaxy stacks. However, 
since the residual "sky" is actually composed of at least 
three different light sources, the background subtraction 
performed in subsection 12.31 has little physical meaning. 
These light sources include "real" sky background, light 
from undetected, and therefore unmasked, galaxies and - 
most importantly - light from the faint outer halo of the 
central LRG. In order to obtain a meaningful background 
measurement we stacked thumbnails at random locations 
within fields from the stacking list, using the same num- 
ber of frames as are in the galaxy stacks and the exact 
same sky level values as determined for individual fields 
in the previous step. For each of the optical bands we 
created one hundred random stacks in every magnitude 
bin to achieve a statistically meaningful distribution of 
the noise characteristics. We then averaged the random 
stacks and subtracted the resulting image from each of 
the magnitude-bin stacks. The background in the LRG 
stack is now defined as all emission in the vicinity of the 
LRG in excess of that of a random, nearby position. 

2.5. Final steps 

As a last step prior to averaging the magnitude bin 
stacks we normalized them with respect to each other 
using the mean apparent magnitude of their input im- 
ages. We then summed the normalized stacks, weighting 
them by their relative total number of stacked frames. 
Finally, for the purpose of absolute photometry calibra- 
tion we matched the final stack to the average appar- 
ent magnitude of the 18.0<mr<18.2 bin, resulting in 
five photometrically calibrated stacks in the u,g,r,i and z 
bands. We note that this method of normalizing magni- 
tude binned stacks is significantly less sensitive to pho- 
tometry errors than normalizing individual images. 

3. STACK ANALYSIS 

Each of the 42579 frames that went into making the 
final stack was obtained with 53.9 seconds of exposure 
time. The final stacks, then, have an effective depth that 
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Figure 4. The empirical LRG group luminosity function (thick 
red dots) was derived by subtracting the luminosity function of ob- 
jects in randomly selected frames (purple dashed line) from that of 
the LRG frames (blue dot-dashed line). The LRG only luminosity 
function is also plotted, along with the best fit power-law curve to 
the faint end slope. 



corresponds to an exposure time of roughly 2.3 Msec, 
equivalent to 40 hours on a 10m class telescope. Obvi- 
ously, such deep images of individual objects would be 
completely dominated by light from neighboring objects, 
flat field errors and scattered light by the PSF. Stack- 
ing a large number of images may currently be the best 
method for studying the faint outskirts of galaxies. 

3.1. PSF stacks 

When analyzing any deep dataset one must be aware 
of the effects of light scatter due to the point spread func- 
tion in the images. Since the dynamic range in the data is 
large, the faint outskirts of the stack can in fact be dom- 
inated by scattered light from the bright galactic center. 
In order to determine the effects of the PSF on the galaxy 
stacks and to properly remove them we stacked synthe- 
sized PSF images created using Robert Lupton's Read 
Atlas Images codeQ. The PSF images were produced 
at the same CCD positions as the stacked LRG images 
for each of the observed bands. However, the synthetic 
PSFs only extend as far as ~12" (58 kpc) in radius and 
do not reach the full extent of the LRG stacks, where 
PSF "co ntamination" may still be importa nt (see, for 
example. Ide Jong|l2008l : lBergvall et al.ll2009f ). We there- 
fore stacked all the bright star images from SDSS in the 
r-magnitude range 8.0< nir <8.2 in each of the five op- 
tical bands. The resulting stacks are saturated inside of 
roughly 5" , where they are unusable, but extend out to 
a radius of 50". We then combined the synthetic pro- 
files with the bright star profiles to create PSF images in 
each band that include the effects of scattered light both 
at small and large radii. The synthetic, bright-star and 
combined PSF i-band profiles are shown in figure [3] along 

^ http:/ /www. sdss.org/dr7/products/images/read_psf.html 



with the stack light profile. In the appendix we perform 
additional testing of the effects of different PSF models 
on our results. 

3.2. PSF deconvolution and Surface brightness profiles 

Several techniques have been proposed and widel y used 
for PSF deconvolut ion from imaging data (e.g. iLucvl 
119741: iHogbomI 119741) . These algorithms, however, typ- 
ically work well either in the inner parts of galaxies or 
their faint outskirts but not both simultaneously. We 
theref ore chose to use t he recently introduced technique 
of Szo moru et al.l ()2010[ ) for PSF deconvolution from our 
LRG stacks. Following this method, we first used Galfit 
3.0 (Peng et a h 2002) to fit a 2D Sersic model with the 
combined PSF images serving as the input convolution 
kernel. For the purpose of pixel weighting, to which the 
GALFIT results are sensitive, we supplied the software 
with the weight masks produced for each stack (subsec- 
tion (531) • We then constructed a synthetic Sersic profile 
from the best-fit parameters to which we added the resid- 
uals that were left after subtracting the PSF convolved 
model from the image. In their paper, iSzomoru et al.l 
(|2010( ) showed that this method is insensitive to varia- 
tions in the Sersic parameters used to create the model. 
In order to increase the flexibility of our model fits we left 
the background level as a free parameter, in effect fitting 
the light profiles with a simple two-component model of 
a constant plus a Sersic function. We also ran GALFIT 
with the background level set to zero in order to test how 
the lack of a constant component affects the model fits. 
Finally we measured the total fiux in the light profile be- 
fore and after deconvolving the PSF from it in order to 
verify that flux is globally conserved. We found that the 
difference between the total flux within the analysis ra- 
dius of 475 kpc is less than 0.5% between the deconvolved 
and the original profiles. 

Figure [5] shows the PSF deconvolved profile of the r- 
band stack (the stack of galaxy images in the r-band) 
along with the GALFIT Sersic model fit. To derive the 
profiles we first applied the IRAF task ellipse on the av- 
erage of the r,i and z stacks while allowing the central 
position, ellipticity and position angle to vary with ra- 
dius. We then used the output table from this fit as 
an input template for obtaining the light profiles of the 
stacks and models. The error bars were derived using 
randomly selected field stacks (see subsection 13.31 for de- 
tails). Also shown in the figure are the initial r-band 
profile and the residual background level as determined 
by GALFIT. 

We note that an ideal deconvolution of the light profile 
would require a radially-varying PSF which is weighted 
by the number of galaxies stacked in each radius bin. 
However, we assume that this would have a negligible 
affect on the properties of the light profiles compared to 
other systematic effects and do not further modify the 
image PSF deconvolution. 

3.3. Profile error analysis 

The unprecedented depth and background uniformity 
of the stacks reveal faint emission hundreds of times dim- 
mer than the typical LRG central surface brightness. In 
images of such low level emission statistical errors are 
not the only significant, and perhaps not even the main. 
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Figure 5. The PSF deconvolved profile of the r-band stack (black 
solid line) is shown along with the best GALFIT Sersic model (red 
dashed line). The sky brightness level as determined by GALFIT is 
also plotted (black dotted line), along with the Icr variation in the 
random stack measurements. The error bars are the standard de- 
viation of a distribution of one hundre d mo del fits using randomly 
selected stacked fields (see subsection 13.31 for details). This pro- 
file is corrected for undetected group galax ies u sing the empirical 
luminosity function described in subsection 13.41 



source of noise. Errors in the faint outskirts of the stacks 
from undetected sources and remaining flattening issues, 
and their affect on the model fits become increasingly im- 
portant but also increasingly hard to asses. We therefore 
used the same random stacks that are described in sub- 
section 12.41 to measure the effects of such biases on the 
stack light profile properties. We subtracted each ran- 
dom background image from its corresponding magni- 
tude bin stack and created one hundred stacks per filter, 
repeating the steps described in subsections 12.31 and 12.51 
We then followed the procedure described in subsection 
13.21 to deconvolve the PSF from the stacks and derived 
a Sersic parameter set for each random stack. Table [1] 
shows that the scatter around the obtained Sersic val- 
ues due to these errors is small, suggesting that the fit is 
weighted toward the luminous inner part of the stacks. 
The Sersic index la deviations of all but the faint u-band 
stack are under 0.1 and the effective radius la deviations 
are under 0.2 kpc. In addition, we used the IRAF task 
ellipse to obtain radial light profiles for each of the ran- 
dom background subtracted stacks, using the same input 
template table as was used in subsection 13.21 The error 
bars in figures [5l [H] and [5] reflect the la scatter of surface 
brightness profile fits to the random stacks. 

Systematic errors in the measured profiles that may 
rise from stacking images of galaxies with a range 
of Sersic parameters see m to be of lesser significance. 
Ivan Dokkum et al] (|2010l . appendix B) stacked hundreds 
of synthetic model images with randomly generated 
Sersic profiles and showed that the stack effective radius 
and Sersic index n match well with the average values of 
the stacked profiles. This suggests that no additional sig- 
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Figure 6. The light profiles of the five stacks are plotted along 
with their best fit model and the residual sky background as deter- 
mined by GALFIT. The quotients of each model and light profile 
pair are also plotted, showing good agreement out to 100 kpc. 
The error bars are the standard deviation of a distribution of one 
hundre d mo del fits using randomly selected stacked fields (see sub- 
section |3]3] for details). 

nificant errors are introduced by our stacking technique. 

3.4. Undetected light correction 

Despite aggressive masking of foreground and back- 
ground objects prior to stacking the data some light from 
undetected, and therefore unmasked, sources remains in 
the images. In order to correct the stack light profiles for 
this effect we used the LRG images to derive an empirical 
luminosity function for the LRG group environment at 
z ^0.34. We started by using SExtractor to extract pho- 
tometry for all the objects in the LRG frames using the 
detection limits that were used for object masking. We 
then binned the resultant values with a bin size of 0.01 
dex in log space and produced a luminosity function (in 
practice we produced an observed brightness function). 
We repeated both steps for an identical number of ran- 
domly selected fields from the same SDSS imaging fields 
and produced a luminosity function for background and 
foreground sources. The LRG group luminosity func- 
tion is then given by the difference between the num- 
ber of sources in a particular brightness bin in the LRG 
fields and the number of random field sources in the same 
brightness bin (see figure |4]). 

The contribution of light from undetected objects to 
the stack profile comes from the faintest group members 
whose luminosity function is expected to be well fitted 
by a power-law function. We therefore divided the LRG 
group luminosity function into an LRG component and 
a power-law component (solid and dotted black curves in 
figure |4]). The best fit power law curve to the data is rel- 
atively shallow, with a slope of $ cx flux~^ '^^ , implying 
that less than 0.3% of the total unmasked light in the 
group (the total light in LRGs and undetected objects) 
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Figure 7. The ellipticity of the rotated stacks as a function of 
radius. The solid curves represent the profiles of the PSF decon- 
volved stacks in the r, i, g and z bands. The ellipticity in all filters 
is relatively fiat out to roughly 60 kpc, where a mild rise is ob- 
served. The initial r-band stack profile (black dashed line) and the 
ILucvl ||T974) deconvolved profile (purple dot-dashed line) are also 
shown. The bottom-right inset shows the r-band rotated stack. 



comes from undetected sources. 

Finally, we extracted a radial light profile from the flux 
of resolved (and masked) objects in the LRG frames and 
fitted it with a single parameter Sersic model. We then 
normalized the best fit curve such that the total flux un- 
der it equaled 0.3% of the total LRG flux. Assuming that 
the radial light distribution of the unresolved sources is 
identical to that of the resolved objects we subtracted 
the resultant correction profile from the stack light pro- 
file. This correction did not change the LRG profile sig- 
nificantly and only moved it within the error bars for any 
given radius bin. 

3.5. Sersic profile fits 

The derived profiles of the g,r,i and z band stacks are 
well fitted by a single Sersic parameter set out to r=100 
kpc, at which radius the profiles deviate from the model 
by 0.2 mag arcsec"^, corresponding to a difference of 
roughly 20%. This is shown in figure [HI where the PSF- 
deconvolved profiles of all five stacks are plotted along 
with their best GALFIT models. The determined Sersic 
index values for the u and g band profiles are n ^4 and 
are slightly shallower than those of the redder colors with 
n ~5 for the r,i and z band stacks. The effective radii 
of all profiles except for the u-band stack are between 11 
and 13 kpc with errors under 0.2 kpc. As might have 
been expected, setting the background level to zero in 
GALFIT increased the best fit effective radii and Sersic 
parameters by up to 10% as GALFIT attempted to fit 
the excess light at large radii. At radii larger than 100 
kpc the profiles deviate from the model and excess light 
is observed in the g,r,i and z band stacks. Figure H] also 
shows that at r > 200 kpc there is 30% to 70% more light 
in the PSF-deconvolved profiles than in their respective 



Table 1 

Stack Sersic parameters 



Filter 


Sersic index 


Effective radius (kpc) 


u 


3.94±f.62 


17.0±9.80 


g 


4.03±0.09 


12.6±0.f9 


r 


5.50±0.05 


13.lit0.f0 


i 


4.86±0.05 


10.9±0.06 


z 


4.9f±0.08 


11.5±0.12 



best fit models. 

The derived profile parameters for the five final stacks 
are presented in tabic [1] 

3.6. Rotated LRG stacks 

Early work in the field of massive galaxy evolu- 
tion showed that the average ellipticity of nearby el- 
lipticals is 0.3 < e < 0.4, with a smaller value for 
the most massive systems (e .g., iSandage et al.l 119701 : 
IBinnev fc de Vaucouleu7slll98lD . The average ellipticity 
of LRGs can be measured from our image stacks. We 
used SExtractor to derive a position angle for each of 
the 42579 galaxies, rotated all frames to a common axis 
and stacked them following the steps described in section 
[21 We then deconvolved the PSF from the stacks using 
the technique discussed in subsection l3.2l Since this tech- 
nique utilizes Sersic parameter fits with radially constant 
ellipticities as the underlying model we teste d our results 
by deconvolving the r-band stack using the lLucvl ([l97l 
algorithm. Figure [7] shows the rotated stack ellipticities 
in four of the bands, excluding the u-band stack due to 
its significantly noisier profile (as can be seen in figure 
El). Also plotted are the original, not PSF deconvolved r- 
band stack (black dashed line) and the Lucy deconvolved 
r-stack profile (purple dot-dashed line) . All profiles were 
obtained using the IRAF task ellipse and all are flat out 
to roughly 60 kpc with an average ellipticity value of 0.26 
between 10 and 60 kpc (e = 0.21 when the centers are 
not excluded). The radial dependence of the ellipticity 
outside of 60 kpc appears to be mfld, with a slight rise 
to e '-^ 0.3 at r = 100 kpc. Furthermore, the proflles of 
both PSF-deconvolved and PSF-deconvolved stacks are 
in good agreement with each other, with less scatter in 
the inner parts of th e former. From a sim ilarly rotated 
stack of BCG images iZibetti et all (|2005| ) found a steep 
rise in profile ellipticity outside of roughly 80 kpc fol- 
lowed by a steep decline at radii greater than 200 kpc. 
In the range 80 to 200 kpc this trend is similar to what 
we observe for LRGs, suggesting that there may be a 
continuum of properties go ing from group to cl uster en- 
vironments. We note that IZibetti et al.l (|2005[ ) did not 
correct their stacks for PSF-induced effects. 

3.7. Color gradients 

It has long been known that the broadband col- 
ors of nearby elli ptical galaxies va r y with radius in 



the optical regime (jVader et al.|[T988l: IPranx et~aLlll989l: 
iPeletier et al.l 11990*^ Nevertheless, studies measuring 
such color gradients typically rely on observations of the 
inner parts of these objects, out to only 2-3 effective radii. 
In figure [51 we plot the color profiles of r-i and g-z. The 
color pairs (r and i, g and z) have similar ellipticity pro- 
file sizes to ensure that the stack depths are compara- 
ble (figure 111) . We compare our results with a study of 
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Figure 8. Broadband color comparison between the inner parts of nearby galaxies from [Peletier et al.l 1)19901 ) and the stacks. The blue 
density plots show the distribution of nearby galaxy colors out to about three effective radii and the red and light blue solid lines follow 
the running means of the sample. The b lue and purple dot-dashed lines represent the r-i and g-z colors of the stacks, respectively, which 
match well with the lPeletier et sH] II1990I ) color profiles. Both colors are flat, within the error bars, outside of about 40 kpc. 



nearby ellipticals bv lPeletier et al.l ()199Cl[ ) and show that 
the profile shapes of the two studies agree well out to a 
radius of ~30 kpc. The stack color values are plotted in 
the ob served frame and are matched to the lPeletier et al.1 
(|1990D rest frame colors by adding a constant. We note 
that outside of about 20 kpc the local sample is composed 
of only a few galaxies where sufficient depth was achiev- 
able. With the advantage of our deep stacks, however, 
we are able to study the colors of LRGs out to more than 
eight effective radii. This is the first time that the col- 
ors of massive galaxies are observed at such large radius, 
providing a new insight on the stars that are found in the 
outskirts of massive galaxies. Indeed, the color profiles 
of LRGs change trend beyond the limit of nearby galaxy 
observations. Despite initially getting bluer in the inner 
40 kpc of both plotted colors, the profiles quickly fiatten 
out to a relatively constant value. The outskirts of LRGs 
are then roughly 0.15 dex and 0.2 dex bluer than their 
centers in the r — i and g — z colors, respectively. 

The color profile is different from that observed by Z05 
at r>20 kpc as Z05 find that BCGs become very red at 
large radii. As we show in the appendix, the Z05 color 
profile may be severely affected by PSF effects at all radii, 
including the stack outer parts. Therefore, the different 
radial color profiles do not necessarily imply that LRGs 
and BCGs are fundamentally different objects. We note 
that the g-r color gradient found by Z05 has a similar 
slope as the g-z profile presented in figure [H The g-r 
profile of Z05 may suffer less from PSF effects than their 
r-i profile. 

4. DISCUSSION 

4.1. Light profiles and the size of massive galaxies 

The first and foremost result that arises from this study 
is that faint, gravitationally bound stellar light can be 
traced in massive elliptical galaxies out to a radius of 
100 kpc. By stacking a large number of faint galaxy im- 



ages we detect light at such distance from the centers of 
massive galaxies with good confidence. In figure [10] we 
show that the total accumulated light at 20< r/kpc <100 
is non- negligible and accounts for roughly 27% of the 
overall flux in the stack. In fact, more than 13% of the 
stack light can be detected at very large radii outside 
of r =100 kpc, or more than 8 effective radii. This is 
especially interesting in light of recent studies that find 
compact massive galaxi es at z ^2, exhibiting effective 
radii of only 1 kpc (e.g..lpaddi et a l.ll2005HTru iillo et al.l 
l2006Hvan Dokkum et al.ll2008D . The size growth of these 
objects is evidently rapid, expanding the physical scale 
of the galaxy by a factor of at least 5 in less than 10 Gyr. 
Unfortunately, the physical growth mechanism cannot be 
directly observed in the LRG stacks as any signal from 
individual galaxies is smoothed and averaged over the 
entire sample. Nevertheless, the lack of a clear change in 
the stack light profile slopes out to 100 kpc suggests that 
the observed light in the outskirts of LRGs is physically 
associated to the galaxies and their inner parts. Fur- 
ther evidence for this comes from the relatively radially- 
independent ellipticity profiles which vary only slightly 
out to 100 kpc. Any other light sources, such as back- 
ground contamination or residual PSF scattering would 
be uncorrelated with the position angle of the LRG as 
measured in individual SDSS images, resulting in a cir- 
cular light distribution. 

Outside of roughly 100 kpc the light profiles of the g,r,i 
and z-band stacks depart from the simple Sersic model 
profile and exhibit excess light (figure [HI). This departure 
from a simple model is observed here for the first time 
in LRGs and it shows that stars at the extreme outskirts 
of massive galaxies follow a different gravitational poten- 
tial than stars in the inner parts. It is known that the 
potential at these radii is dominated by the properties of 
the dark matter halo, implying that the light profile is 
not necessarily expected to follow the same model that 



8 



Tal & van Dokkum 



20 - 



22 - 



CT" 24 
o 

in 

ig 26 
c 



CT> 



(J 



28 - 



30 - 



32 



r-bond stock profile 
n = 5.5 model fit to stock 
ICL+BCG profile (Zibetti 05) 
n=4 fit to ICL+BCG profile - 



v.. 



\\ < 



10 



Radius [kpc] 



100 



1000 



Figure 9. A comparison between the light profiles of our LRG r- 
band stack and the ICL profile from : 2ibettl~ejral.| ^005). The ICL 
profile departs from a single parameter Sersic model at 50 kpc, or 
double the departure radius of 100 kpc that is observed in the LRG 
stack. This suggests a more significant population of intorgalactic 
stars in massive clusters than in groups. 



describes the inner stellar body. Alternatively, this ex- 
cess light may simply be the residual background in the 
images, reflecting unresolved light from the group envi- 
ronment in which LRGs typically reside. 

Excess light was also observed by Z05, who studied 
the ICL around brightest cluster galaxies from a stack of 
683 SDSS images. Such galaxies typically live in dense 
halos with total mass of 10^"* to IQ^^Mq and are inher- 
ently different from LRGs whose group halos are a few 
times lO^^Af0 in mass. Z05 found that in clusters, this 
"extra light" constitutes only a small fraction of the to- 
tal cluster profile, accounting for less than 11% of the 
light inside of 500 kpc. Nevertheless, the ICL profile de- 
parts from a single parameter Sersic model already at 
r ~ bOkpc, compared to the departure radius of 100 kpc 
that is observed in our LRG stacks (figure [9]). This sug- 
gests that the massive clusters studied by Z05 may more 
readily support a population of intergalactic stars than 
the groups in which LRGs reside. In their paper Z05 
correct their light profiles for unresolve d cluster sources 
using the luminosity function given by iMobasher et al.l 
([2003). We note that the PSF, which is not deconvolved 
from the ICL+BCG profiles presented in Z05 may scat- 
ter light at all radii and increase the errors of the Sersic 
model fit. 

Unlike the outer parts of LRGs, the centers of these 
galaxies are not well resolved in our stacks. Studies uti- 
lizing high resolution HST images showed that the pro- 
file at the inner parts of nearby ellipticals often departs 
from the Sersic model that traces their outskirts. More 
specifically, the most mass ive ellipticals exhibit fiattened 
central light profiles (e.; 
1991 'Lauer et al.' 19951; 
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Kormendv et al., (,2009i) used a compilation of HST and 



Figure 10. The radial stellar light fraction of the r-band stack 
in four radius bins of 20,100,200 and 300 kpc. Note that a non- 
negligible fraction of the light is detected in the extreme outskirts 
of the stack. 



ground based data to show that although well fitted by a 
Sersic model out to large radii, the most massive Virgo el- 
lipticals exhibit 1 kpc scale cores. In our stacks we cannot 
resolve such physical scales as 1 pixel in the SDSS data is 
equivalent to 1.9 kpc at the stack mean redshift of 0.34. 
We are nevertheless able to confirm the excellent fit of 
massive elliptical galaxy profiles to a single Sersic profile 
out to a few effective radii that IKormendv et al.l ()2009l ) 
found for individual Virgo galaxies (reaching A/x^ > 0.2 
mag arcsec^ at r\ > 100 kpc). 

4.2. How much light is missed? 

The deep stacks allow us to test how much light is 
missed in typical studies of the profiles of individual 
LRGs and derive a correction factor that can be applied 
in such cases. To do so we first selected all the LRGs in a 
single magnitude bin, 18.0<mr<18.2, and used GALFIT 
to produce a Sersic model to each object individually. We 
then excluded all fits with errors of more than 10% in ei- 
ther the n parameter or the effective radius, resulting in 
a mean effective radius value of 11.7 kpc. The difference 
between this value and the one derived by GALFIT for 
the stacked image {r^ =13.1 kpc) is then ~10%. This 
implies that surveys may underestimate the size of mas- 
sive red galaxies by this amount. The total flux in the 
stack, however, accounts for ^20% more than the mean 
value for the individually derived profiles, suggesting that 
a non-negligible amount of light is typically missed and 
that the total stellar mass is underestimated. 

4.3. Minor mergers and the LRG color profile 

It has long been known that the color profiles of 
nearby massive ellipticals exhibit a relatively smooth gra- 
dient toward bluer colors from the galaxy centers out- 
ward. Line in dex measurements fe.g..l Carollo et a"Llll993l : 
iDavies et al.l [1993; .Spolaor et al.. .20101 ) and studies of 
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high redshift galaxies (jTamura et al.|[2000l ) provide evi- 
dence that this observed color gradient originates from a 
radial slope in stellar metallicity. It has also been claimed 
in these studies that age plays only a secondary role in 
producing this trend. Nevertheless, it is yet to be de- 
termined which physical mechanisms create and evolve 
the observed color gradient. With the color profile in 
mind we will discuss in this section a previously suggested 
model in which minor mergers and low mass accretion 
events are the main supplier of stars to the outskirts of 
massive galaxies. We will also utilize our LRG stacks to 
provide additional pieces of evidence that support this 
model. 

4.3.1. Creating the faint stellar halo 

It has been suggested that the size and mass growth of 
massive galaxies are dominated by minor mergers which 
contribute a relatively constant flow of accreted mass 
over a lon g time. This is evideii t from analytic calcula- 
tions (e.g. JBezanson et al.ll2009D. numerical simulations 
of lar ge volumes (e.g.. lLucia fc Blaizotll2007t iNaab et aD 
2007f) and observations of nearby e lliptical galaxies (e.g., 
van Dokkun]||2005HTal et al.ll2009[ ). The contribution of 
accreted mass, however, is not expected to be uniform 
throughout the galaxy. This can be shown using very 
simple energy arguments following the assumption that 
accreted mass is more likely to stay near the radius at 
which is was accreted if its total energy budget is sim- 
ilar to the gravitational energy of the accreting galaxy 
at that radius. In other words, if the escape velocity of 
a star in an infalling galaxy equals the orbital velocity 
of the accreting galaxy at the radius at which the star 
escapes, this star will likely fall into orbits close to that 
radius. If on the other hand the star escapes the in- 
falling galaxy at any other radius it will either fall closer 
to or be thrown farther from the center of the accreting 
galaxy, thus contributing to a more uniform mass distri- 
bution. This argument can be described by the following 
equation: 

KTc (1) 



90 



M 
aire 



where V^^.^ is the circular velocity of the accreting galaxy 
and Vg™^ is the escape velocity from the infalling galaxy. 
This can be expressed in terms of the masses of the ac- 
creting and infalling galaxies, M and m and the accretion 
and escape radii, Ra and r, 



Ra — 



M r 
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Equation [5] implies that more massive infalling galaxies 
will preferentially deploy most of their stars closer to the 
center of the accreting galaxy. However, the mass ratio 
M/m cannot be close to 1 since such major merger will 
violently disrupt both galaxies and not simply disperse 
the stars of one into the stellar body of the other. In 
addition, if the mass ratio is high, the contribution of 
the infalling galaxy to the colors of the accreting galaxy 
will be minor. We therefore assume typical values of 
4<M/m<10 for the mass ratio and 2< r/kpc<10 for the 
escape radius, both roughly correspond to the accretion 
of an L* type galaxy. It is therefore convenient to re- 
write equation [2] using the suggested average values for 
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Figure 11. A comparison between the position angle of stellar 
light a nd that of ti dal feature light in galaxies from the OBEY 
survey I ITal et al.l2d09 ') . The correlation between the orientation of 
tidal features and the stellar body suggests that minor interactions 
cannot be ruled out as a source of stellar light to the outskirts of 
massive galaxies. 
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This simple model provides a possible explanation for 
the observed size growth from z=2, as well as the smooth 
gradient and flattening in the stack color profile. It en- 
ables infalling low mass galaxies to easily deploy the ma- 
jority of their stars outside of the accretion radius Ra and 
it suggests that scattered stars that are stripped from 
their galaxy at any other radius may end up close to the 
center. This may be the case if the accreting galaxy is 
indeed red and if a significant part of all infalling galax- 
ies are bluer than the LRG. In addition, this scenario 
supports minor mergers with large M/m values as the 
main size growth mechanism as stars are preferentially 
being deployed far from the center of the accreting galaxy 
without increasing its total mass significantly. 

4.3.2. The frequency of minor mergers 

The frequency of low mass accretion events is a critical 
factor in assessing the importance of these interactions 
to the size and mass growth of massive ellipticals. Minor 
mergers must not be rare in order for the compact galaxy 
at z=2 to experience the observed rapid size growth. We 
now also know that the light at r > 20 kpc comprises 
roughly 40% of the total stellar light (figure [TU]), sug- 
gesting that a constant stream of accreted mass is re- 
quire d. Th is implied mass growth confirms the results 
from Ivan Dokkum ct al.. (,2010i) who found that nearby 
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ellipticals are roughly twice as massive as compact z—2 
galaxies. 

Observations of early type galaxies reveal ample evi- 
dence for close gravitational interactions between ellip- 
ticals and their neighboring galaxies. Such interactions, 
mainly in the form of minor mergers, leave their signa- 
ture on the stellar bodies of the massive galaxy for a time 
period that depends on the interaction itself and that can 
last for up to a few Gyrs. Initial attempts to characterize 
tidal features around nearby ellipticals relied heavily on 
qualitative descriptions of the remnant morphology but 
neverthe less revealed evidence for interaction in many 
systems (jSchweizer fc Seitzeilll992|) . More recent studies 
quantified an interaction parameter by measuring excess 
light in the tidal features compared to smooth models of 
the systems. Such studies analyzed volume limited sam- 
ples of nearby galaxies, finding that more than 70% of 
all ellip ticals have had one or more recent m inor merger 
events ()van DokkumI I2OO5I : iTal et al.l 120091) . The mass 
growth that is inferred from this high interactio n rate can 
be es timated at a factor of 2-3 since z=2 (see I Tal et all 
[2009l for details). 

4.3.3. Minor mergers and the ellipticity profile 

Assuming that the distribution of infalling galaxies 
and therefore also tidal features is spherically symmetric 
around the elliptical, the stack outskirts are expected to 
be relatively round, in apparent confiict with the trend 
shown in figure [71 It is however possible that gravita- 
tional interactions contribute to the ellipticity of the host 
galaxy or alternatively occur preferentially along the ma- 
jor axis of the dark matter halo, thus aligning the resul- 
tant tidal features with the position angle of the ellip- 
tical. In order to test this we examined galaxies from 
the OBEY survey (Tal et al. 2009), a complete sample 
of nearby ellipticals that was utilized to study and quan- 
tify tidal features in the stellar bodies of these objects. 
We measured the overall ellipticity of each galaxy with 
a tidal parameter value of at least 0.07 and compared it 
to the second moment of light distribution in the model- 
subtracted residual image. By measuring the second mo- 
ment of residual light distribution we quantified, in effect, 
the position angle of tidal features. Figure [TT] shows that 
a correlation indeed exists between the orientation of the 
stellar body and that of the tidal features. A Spearman's 
rank test finds that the probability that this correlation 
was drawn from a random distribution is less than 0.1% 
We also note that only two galaxies exhibit /S.PA > 60° . 
Although such a suggestive correlation does not strongly 
support any single scenario as the main mechanism for 
creating the observed faint halos, its existence hints that 
minor mergers at least cannot be ruled out as a signif- 
icant contributor to the stellar bodies of LRGs at large 
radii. 

This analysis is therefore consistent with the idea that 
minor mergers and gravitational interactions likely play 
an important role in determining the properties of mas- 
sive red galaxies. Foremost is the blue color index of the 
profile at r > 40 kpc which suggests that the outer halo 
is composed of younger or alternatively more metal poor 
stars compared to the center. This probably means that 
the stellar populations of the outskirts were formed sep- 
arately from those in the center and probably accreted 
at a later time. This scenario is supported by the ob- 



served high rate of tidal features around nearby elliptical 
galaxies, although we note that the LRGs are typically 
more mas s ive by a factor of >2 than the galaxies in the 
ITal et all (|2009l ) sample. Our analysis shows that such 
accretion events likely deploy most of the stars at large 
radii. 

5. SUMMARY AND CONCLUSIONS 

In this paper we stacked more than 42000 images of 
luminous red galaxies in order to study the faint light of 
these objects at large radii. In our stacks we detected 
stellar light out to radii greater than 100 kpc, thus pro- 
viding a correction factor to the true size and overall 
stellar mass of LRGs. This is the first time that such 
an analysis is performed using a dataset of this scale, 
reaching unprecedented depth for z > 0.01 galaxies. The 
relatively flat ellipticity profile (figure [7]) verifies that the 
light detected at large radii is physically associated with 
the stellar body. Interestingly, the profiles suggest in- 
creased ellipticity at large radii of the average LRG pro- 
file. 

In agreement with iKormendv et al.l (|2009[ ). we con- 
firmed that on average, the light profiles of massive el- 
lipticals can be well described by a single Sersic model 
out to roughly eight effective radii. Outside of 100 kpc 
the profiles deviate from a simple Sersic model and ex- 
hibit extra light in the r, i and z-band stacks. This excess 
light can probably be attributed to unresolved intragroup 
light or a change in the light profile itself. Differentiat- 
ing between these possibilities, however, may be difficult 
as both can have a similar effect on the profile shape at 
large radii. 

Finally, we utilized the five optical bands of SDSS to 
study the colors of these galaxies and showed that the 
well known decrease in color index out to 2-3 effective 
radii flattens out and stays blue compared to the galac- 
tic centers out to the detected stack limit. Although this 
finding by itself does not favor any one stellar population 
evolution scenario, it suggests that the central 20 kpc 
evolve somewhat differently from the rest of the galaxy. 
Previous studies of line indices in early type galaxies sug- 
gest that this difference can probably be attributed to a 
difference in metallicity. 
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APPENDIX 
THE WINGS OF THE PSF PROFILE 

Light scattered by the wings o f the PSF profil e was shown to contribute noticeably to the color profiles of the outskirts 
of spiral galaxies. In his paper, Ide Jond ()2008[ ) argues that the choice of PSF model used for image deconvolution is 
crucial for analyzing data of spiral galaxies. The author further claims that the red stellar halos found by previous 
studies were artificially produced by PSF models that were poorly sampled at large radii. In this appendix we confirm 
this observation and show that the widely used synthetic PSF images may lead to inaccurate color gradients even in 
the peaky light profiles of massive red galaxies. 

In order to test the effects of PSF profile selection on the properties of the derived Sersic model we followed the 
procedure described in subsection 13.21 to deconvolve the stacks using two different PSF images in each observed band. 
The first profile is identical to the one used throughout this study and it combines a synthetically produced PSF model 
with a bright star stack (subsection 13. 1|) . while the second utilizes only the synthetic PSF image produced by the Read 
Atlas Image code. While the u,g,r and z models are less than 15% larger when using only the synthetic PSF image, 
the i-band stack exhibits an increase of more than 25% in this case. The increased growth in model size in the i-band 
compared to the other stacks implies that the halo of these galaxies may erroneously appear to be red. This is most 
clearly evident in figure [T^ where we plot the r-i color profiles of the two PSF deconvolved models and show that 
in the synthet ic only case an artificial red halo appears outside of roughly 40 kpc. We conclude, in agreement with 
Ide Jond ()2008f ) . that a proper choice of well sampled PSF model is critical for studying the faint outskirts of galaxies. 
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Figure 12. Color profile comparison between the combined PSF deconvolved stack (blue line) and the synthetic only PSF deconvolved 
stack (red line). The shaded areas show the la error bars derived using randomly selected field stacks (subsection 13.31 1. The use of 
insufficiently sampled PSF profile artificially creates an red halo in our stacks. 



